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For odd anharmonic oscillators, it is well known that complex scaling can be used to determine 
resonance energy eigenvalues and the corresponding eigenvectors in complex rotated space. We 
briefly review and discuss various methods for the numerical determination of such eigenvalues, 
and also discuss the connection to the case of purely imaginary coupling, which is PT-symmetric. 
Moreover, we show that a suitable generalization of the complex scaling method leads to an algorithm 
for the time propagation of wave packets in potentials which give rise to unstable resonances. This 
leads to a certain unification of the structure and the dynamics. Our time propagation results 
agree with known quantum dynamics solvers and allow for a natural incorporation of structural 
perturbations (e.g., due to dissipative processes) into the quantum dynamics. 
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I. INTRODUCTION 

Quantum tunneling is one of the intriguing phenomena in quantum mechanics. Since the seminal work of G. Gamow 
O r on the theory of the alpha decay of a nucleus [1], quantum tunneling of a particle trapped in a metastable potential 
well has been studied in detail in many areas of physics. Indeed, in mathematical terms, the decay width of the 
corresponding resonance state in an, e.g., cubic potential can be traced back to an instanton configuration which is 
j3 ■ a solution of the classical equations of motion of a particle moving in the "inverted" potential —V(q) in such a way 
that the classical action of the particle along its trajectory remains finite even though the time domain covered by 
its trajectory is the entire space R. Here, the instanton configuration covers the domain naturally associated with 
the tunneling process from a relative minimum of the potential to a point where a horizontal line would emerge from 
the "tunnel" potential. An approximation to the width can be obtained by considering the fluctuations around the 
' classical path, and by a subsequent evaluation of the Fredholm determinant describing the fluctuations. This is quite 
(~^. \ analogous to the case of double- well-like potentials [2, 3], where the instanton configuration describes an oscillatory 
■ motion covering a trajectory oscillating between the two degenerate maxima of the inverted potential. 

At the same time, the structure of the potentials depends very much on the complex phase of the coupling constant. 
Let us assume a metastable potential approximated by a one-dimensional harmonic oscillator perturbed by a cubic 
term yfg q 3 , where g is the coupling constant. For purely imaginary y/g = i/3, with (3 € R, the Hamiltonian is invariant 
t^" under the composed application of a parity transformation and a time reversal and thus called 'PT-symmetric. In 
some sense, the PT-symmetric case is a natural generalization of hermiticity (or even self-adjointness) to the complex 
, domain, and it will be verified here that a number of concepts known to be applicable to "stable" anharmonic 
oscillators are applicable to the PT-symmetric cubic potential as well. 

Coming back to the case of real coupling parameter ^/g, we note that known methods for the numerical determi- 
nation of resonance energies and widths include the diagonalization of a complex rotated Hamiltonian matrix, Borel 
resummation [4-6] in complex directions of the parameters, and strong-coupling expansions. These three methods 
will be discussed and contrasted here as they are important for the construction of an adiabatic, complex transformed 
time propagation algorithm that is also described in this article. Namely, we attempt to solve the problem of how to 
propagate a wave packet that moves under the influence of a potential with metastable resonances and which may 
thus even "escape" to the classical region of attraction where the potential assumes large negative values, without 
the need for a temporally adjusted numerical grid and without the need for the introduction of transparent boundary 
conditions. Moreover, we attempt to construct this algorithm using complex resonance eigenstates, thus giving a 
manifest interpretation to the complex resonance state and energies within the time propagation method, including 
the back-transformation of the complex rotated and propagated states to the normal coordinate representation. 

This paper is organized as follows. In Sec. II, we recall some basic definitions related to the cubic potential, to 
perturbative and strong-coupling expansions and to the method of complex scaling. Also, corresponding numerical 
investigations are discussed. In Sec. Ill, we discuss a method for time propagation in the cubic potential, which leads 
to the above mentioned desired unification. Finally, a summary is presented in Sec. IV. An appendix is devoted to 
the discussion of potential applications of the algorithms discussed here within a solid-state physics context. In the 
entire article, we attempt to follow a rather detailed style in the presentation and hope that the reader will not find 
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FIG. 1: (color online.) Energy dependence of the function B(E,g) for the VT- 
symmetric case. Results are presented for the coupling parameters ^fg = i/3 with 
P = 1/10 (red circles), /3 = 1/5 (green squares) and /3 = 1 (blue triangles). 



the level of detail excessive. 



II. COMPLEX RESONANCE ENERGIES VERSUS REAL PT— SYMMETRIC SPECTRUM FOR THE 

CUBIC OSCILLATOR 

A. Orientation 

Although our analysis can easily be generalized to quintic and other "odd" potentials with metastable resonances, 
we start for simplicity with the one dimensional Hamiltonian of the cubic potential in the normalization 

which recovers the harmonic oscillator spectrum Em = in the limit of a vanishing coupling constant g — > 0, where 

N is a nonnegative integer and constitutes the quantum number of the state. For nonvanishing g, the operator (1) 
is not self-adjoint, and it does not possess a spectrum of discrete, real energy eigenvalues. For real and positive g, 
the cubic Hamiltonian (1) possesses resonances, i.e. poles of the resolvent (H c — E)^ 1 which is a priori defined as a 
holomorphic function for ImE > 0. The poles (resonances) become apparent in the meromorphic analytic continuation 
of the resolvent to ImE < 0. The resonance energy eigenvalues are of the form 

E N (g)=Re[E N (g)]-i^l. (2) 

By contrast, for purely imaginary coupling ^fg = i/3, the Hamiltonian (1) is 'PT-symmetric, i.e. it is invariant under 
a simultaneous parity transformation q — > —q and a time reversal operation t — > —t, the latter being equivalent to an 
explicit complex conjugation of the Hamiltonian and thus to the replacement (3 — > —(3. 

Based on numerical evidence, it has been conjectured around 1985 by D. Bcssis and one of us (J. Z.-J.) in private 
communications that the spectrum of the PT-symmetric Hamiltonians of odd anharmonic oscillators should consist 
of real eigenvalues, even if these Hamiltonians are obviously not self-adjoint. C. M. Bender and others have recently 
studied PT-symmetric Hamiltonians quite intensively (see, e.g., Refs. [7-11]). Note, in particular, that even the 
quartic anharmonic oscillator with negative coupling becomes a PT-symmetric Hamiltonian with real spectrum when 
endowed with appropriate boundary conditions imposed on the wave function as a function of a complex coordinate, 
as detailed in Eq. (4) of Rcf. [12]. The conjecture has recently been supported on mathematical grounds (see Refs. [13- 
16]). Important further contributions to the development of the theory of 'PT-symmctric quantum mechanics have 
been summarized in Refs. [9-11]. 
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TABLE I: Ground-state energy of the cubic Hamiltonian for the 'PT-symmetric 
case y/g = i /3. Results have been obtained by the diagonalization of the complex- 
transformed Hamiltonian in the basis of the harmonic oscillator wavefunctions ( "ex- 
act" values) and by solving the quantization condition given in Eq. (7). The appar- 
ently "converged" decimals for the solutions of B(E,g) = | are underlined. 






E (g) 


Solution of B(E, g) = \ 


1/10 


0.512 538 145 


0.512 538 145 


1/8 


0.518 760 345 


0.518 760 344(1) 


1/6 


0.530 781 759 


0.530 781 77(1) 


1/4 


0.558 372 124 


0.558 377(4) 


1/2 


0.645 877 080 


0.645 7(3) 



B. Imaginary coupling parameter and real energies 

The discussion in this section is a priori relevant only for the case yfg = i/3, with f3 G R, but we will keep a general 
g in the first considerations, because the general formulas (as a function of g) will be useful later in this article. First, 
we recall that, as explained in Refs. [2, 3], the spectrum of anharmonic oscillators can be described in many cases 
by two functions B and A. Respectively, these are related to the perturbative expansion about the minima of the 
potential and to the tunneling of the quantal particle from one degenerate or quasi-degenerate minimum (if it exists) 
to the other degenerate or quasi-degencrate minimum of the potential. Specifically, around Eq. (3.31) of Sec. 3 of 
Ref. [2], it is explained how, in the case of intcgrable systems, a "perturbative B function" can be obtained from the 
perturbative expansion of the logarithmic derivative of the wave function, which fulfills a differential equation of the 
Riccati type and which is subject to a uniqueness condition which gives rise to the integer N on the right-hand side of 
Eq. (4). Using techniques outlined in Ref. [17], the function B(E,g) can be easily determined for the cubic potential, 
and the first terms read as follows, 

/ 7 15 2 \ o /1365 1155 ,\ m , , N 
B(E,g) = E + g(- + -E^ + g^—E + —E^ + 0(g*). (3) 

By formulating the initial perturbation as y/gq 3 , we have obtained a perturbative expansion for B which involves 
integer powers in the coupling constant. As usual, the perturbative expansion for the A^tli energy level can be obtained 
by inverting the condition 

B(E,g) = N+±, (4) 

and we obtain the standard Raylcigh Schodingcr perturbation theory (RSPT) series for the A^th level of the cubic 
potential, 

oo 

E N (g) = J2^N,K9 k , (5) 
k=a 
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where K is the order of perturbation theory, and the leading pcrturbative coefficients read as follows, for a general 
level with quantum number TV, 
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The minus signs are explicitly indicated to illustrate that all pcrturbative coefficients are negative for all levels, except 
for the leading term N + ^, which stems from the unperturbed harmonic oscillator. 

If perturbation theory determines the energy eigenvalues completely in the 'PT-symmctric case and there are no 
instanton configurations to consider, then the quantization condition (4) can be reformulated as 
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B(E,g)) 



(J. 



(7) 



This quantization condition is formulated such as to display the analogy to those for more complex potentials like the 
double well [see Eq. (2) of Ref. [18]], but it lacks the "instanton A function" present in the cited equation. A Bohr- 
Sommerfeld quantization condition of the form (7) is relevant for stable anharmonic oscillators like, e.g., the quartic 
one with a perturbation proportional to gq 4 for positive coupling g, where there are no instanton configurations 
to consider and therefore no "instanton A function" present. We now investigate to which extent this quantization 
condition could be relevant for the cubic anharmonic oscillator. 

To this end, we first recall that the cubic Hamiltonian (1) displays PT-symmetry for imaginary coupling y/g = i/3. 
The spectrum of the cubic oscillator becomes real in this case, as a consequence of pseudo-Hermiticity [19], and the 
perturbation series (5) becomes an alternating, factorially divergent series in the variable g = — /3 2 ; series of this 
type are typically Borcl-summablc. Wc are therefore led to the conjecture that the quantization condition (7) should 
describe the energy levels of the cubic Hamiltonian (1), but only for the PT-symmetric case ^fg = i/3. 

As discussed in Sec. II and as shown in our previous paper (Ref. [20]), the direct analysis of the perturbative B 
function which enters the quantization condition (7), without any detour via the perturbation series (5), can be used 
to investigate our conjecture formulated above, in order to calculate the energy eigenvalues of the cubic Hamiltonian 
(1) directly. Here, we explore this approach in the case of imaginary coupling ^fg = i/3. To this end, we interpret 
B(E,g) as a function of E (at fixed g) and numerically determine the points at which this function assumes values 
of the form N + ^. According to Eq. (7), these points correspond to the real energy eigenvalues Epj(g) of the cubic 
Hamiltonian in the PT-symmctric case. In many cases, direct investigation of the pcrturbative B function has proved 
to be numerically more stable than that of the Borel-resummed series (5) (see alse Ref. [20]). 

In Fig. 1, for example, we display B(E,g) as a function of the (real) energy argument E for different fixed values 
of the coupling parameters g. Specifically, we consider the cases y/g = i/3 = i/10 (circles), ^fg = i/3 = i/5 (squares) 
and yfg = i/3 = i (triangles). The points where B(E,g) = N + \ are clearly displayed. The (real) ground-state 
energy E^{g) for the PT-symmetric case is presented in Table I and compared to reference data obtained via a 
diagonalization of the Hamiltonian matrix. The ground-state energy is well reproduced at /3 = 1/10 (up to 9 decimal 
digits). For stronger coupling, there is a larger uncertainty in the determination of the ground-state energy value 
T/v = o(g) because the power series (3) is divergent for all nonvanishing g, and, hence, resummation techniques are 
required for its calculation. In all cases, the function B(E,g) is calculated by means of a generalized Borel-Pade 
method, similar to Ref. [20]. We observe, in higher (Borel) transformation orders, an oscillatory behaviour of the 
Borcl integral, evaluated along complex directions, for larger g. These oscillations cannot be overcome when the 
transformation order is increased and represent a fundamental limit of the convergence of resummed weak-coupling 
perturbation theory in the case of a large (modulus of the) coupling parameter g. All numerical experiments support 
the conjecture (7) for the PT-symmetric case. 



C. Real coupling parameter and complex resonance energies 



We again consider the Hamiltonian (1), but this time for real and positive g. First, let us note that the structure 
of the perturbation series (5) is of course unaffected by a change in the complex phase of the coupling parameter. 



5 



0.15 



0.1 



0.05 - 



X 



o.o 



-0.05- 



-0.1 



-0.15 





t 

1 

I s>*(xd) 1 

\ / 

V U( X ) j 

V / 
\. / 
^\ / 


/ 

/ 

/' 

/ 

t 





-0.3-0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

X 

FIG. 2: Instanton configuration for the cubic potential. The dashed line is the poten- 
tial U(x) = X 3 ~~ \ X 2 - The solid line is the worldline of the instanton configuration 
Xd(i) — [cosh(f) + l] -1 , which reads in inverted form t(xd) = ±arccosh[(l — Xci)/Xd]- 



However, we now have instanton configurations to consider. Let us consider for a moment the classical Euclidean 
action, corresponding to (1), 



%(<)] = / dt 



1 fd 



1 



2 \dt q{t) I +n9( t r+V9q(t) 



(8) 



This action describes the motion of a particle in the potential ~^/gq(t) 3 — \ q(t) 2 . Via the change of variable 
q(t) = —x(t)/y/9i we obtain the action, 



%(*)] = - / dt 



1 fd 



2 [g- t x(t)) +^x(ty-x(ty 



(9) 



for which the (redefined) potential now reads U(x) — X 3 ~ h X 2 - The (classical) instanton configuration is (see Fig. 2) 



Xd(t) 



cosh(t) + 1 



5[ X ci(t)] = ^ , 



(10) 



which fulfills U(xd(t = 0)) = f/(|) = 0. An integration about the fluctuations around the instanton path in the 
partition function using methods described in Ref. [21] then leads to the known result 



/7T g 



exp 



2 
155 



(11) 



for the imaginary part of the ground-state energy. Observe that in contrast to the quartic potential, where two 
degenerate instanton configurations are present because of the reflection symmetry of the quartic potential [21], the 
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cubic potential has no reflection symmetry and only one instanton. We conclude that the "instanton A function" for 
the cubic oscillator should have the leading terms 

A(E,g) = ^-+0(g). (12) 
lbg 

However, we stress that the analogue of the quantization condition (7), suitably generalized for odd anharmonic 
potentials, has not yet appeared in the literature to the best of our knowledge. This problem is currently under 
investigation. It is clear that a suitable quantization condition should involve the A function in such a way that the 
decay rate follows naturally by an expansion in both analytic and nonanalytic terms for g small, so that a so-called 
resurgent expansion [22-24] is obtained which allows for the nonanalytic behaviour of the imaginary part of the 
resonance energy as g — ► + . 

If the (generalized) quantization condition for g > were known, then we could use it in order to calculate complex 
resonance energies, via a direct resummation of both the pcrturbative B and the instanton A function, in a similar 
way as was done for the 'PT-symmetric case in the previous section. We recall that for the PT-symmetric cubic 
potential, where we resummed only the B function, and for the Fokker-Planck potential (see Ref. [20]), where we 
resummed both the B as well as the A function, the direct resummation of the quantization condition did not give 
a satisfactory numerical accuracy for the energy levels for moderate and large values of the coupling constant. Here, 
our goal is to describe a numerical method which is applicable to all domains of the coupling constant, and we thus 
continue here with a comparison of three methods for the calculation of resonance energies, including an evaluation 
of the domains of their respective applicability. 

Method I. We apply complex scaling (see Refs. [25, 26]), q — ► ge to the cubic oscillator, which results in the 
Hamiltonian 

H c (0) = c-™ + I ^ + V?9 3 c 5i9 ) . (13) 

The diag^onalization of this complex scaled operator is carried out in the basis of harmonic oscillator wavefunctions 
{ ( t > n(,q)} n =Q l , for large enough iV max , and the variation of the resonance energies under a suitable increase of iV max is 
used to investigate the numerical uncertainty of the results. This allows us to numerically determine the (complex) 
resonance energies of the original cubic Hamiltonian (1). As discussed in Refs. [26, 27], these resonance energies 
are independent of 8, provided we choose 6 sufficiently large so that the rotated branch passes the position of the 
resonance under investigation. 

Method II. We resum Rayleigh-Schrodinger perturbation theory (RSPT) in complex directions of the parameters. 
To this end, we use the standard RSPT series for the cubic potential as given in Eq. (5). Notice that the RSPT 
series is nonaltcrnating in integer powers of g for real g. We then employ the Borel-Pade summation method with 
the (Laplace) integration in the complex plane, as given by Eqs. (196)-(198) of Ref. [28]. The method has also 
been discussed in detail elsewhere [28-30], and it has been put on rigorous mathematical grounds recently [31] in the 
framework of distributional Borcl summability. We employ an integration contour (called C+i in the conventions of 
Ref. [28]) which leads to a negative sign of the imaginary part of the resonance energy eigenvalue, consistent with 
Eq. (2). The accuracy obtainable using this weak-coupling method is restricted by oscillations of the transforms in 
higher orders, which are analogous to those reported in Ref. [20] for other applications of the Borel-Pade transforms. 
Indeed, at a relatively moderate coupling g = 0.6, the ground-state energy as determined by resummation cannot be 
calculated to better accuracy than 

E Q {g = 0.6) = 0.554(1) - 0.351(6) i (14) 

by resummation. Note that the oscillations cannot be overcome when the (Borel) transformation order is increased 
and represent a fundamental limit of the convergence of resummed weak-coupling perturbation theory in the case of 
a large (modulus of the) coupling parameter g. A more accurate result obtained by method I for the same coupling 
parameter is 

E (g = 0.6) = 0.554 053 519- 0.351 401 778 i. (15) 

Method III. We employ a strong-coupling expansion in complex coordinates, to complement the weak-coupling 
perturbative method. Specifically, we employ the so-called Symanzik scaling q — > qg^ 1 ^ 10 and rewrite the cubic 
Hamiltonian into a scaled one H c — > H s with the same eigenvalues but a fundamentally different structure, 



H s = g^ (H e + £g- 2 /^ 
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TABLE II: The coefficients Ln,k of the strong-coupling expansion (17) of the eigenvalues Ei^(g) of the cubic 
potential. For K = 1, the coefficients are zero. For K — 3, based on numerical evidence, we conjecture that the 
coefficient is independent of N, purely real and equal to 1/108 = 0.009259259 . . . 



K_ 7V = 

0.617 160 050 - 0.448 393 023 i 

2 -0.013 228 193 + 0.040 712 191 i 

3 0.009 259 259 + 0.000 000 000 i 

4 -0.000 294 361 - 0.000 905 951 i 



N = 1 

2.193 309 731 - 1.593 532 797i 

-0.022 015 998 + 0.067 758 274 i 
0.009 259 259 + 0.000 000 000 i 

-0.000 141 177 - 0.000 434 499 i 



N — 2 

4.036 380 020 - 2.932 601 744 i 
-0.027 024 360 + 0.083 172 425 i 
0.009 259 259 + 0.000 000 000 i 
-0.000 118 189 - 0.000 363 747 i 




Coupling parameter g Coupling parameter g 



FIG. 3: (color online.) The real (left panel) and the imaginary part 
(right panel) of the exact resonance energy of the ground state for the 
cubic potential are displayed as a function of g (solid lines). The real 
part approaches the values | as g — > 0, whereas the imaginary part 
vanishes in this limit. The exact numerical values are compared to the 
sum of the first four leading terms of the strong-coupling asymptotics 
(dashed lines) for large coupling parameter g. The strong-coupling 
asymptotics are defined in Eq. (17), and the coefficients are listed in 
Table II. The first four terms of the strong-coupling expansion approx- 
imate the exact resonance energies up to surprisingly small values of g, 
but both the real as well as the imaginary part deviate substantially 
for g < 0.025. 



A strong-coupling perturbation expansion can thus be written for each energy Ejsr(g), which reads 

oo 

En{ 9 ) =9 1/5 J^L N , K g- 2K / 5 , (17) 

where Ljv,o is just equal to the Nth level of Hi. Based on the strong-coupling expansion, it is easy to see that a 
rotation angle 6 = ir/5 = 36° is sufficient to uncover all resonance energies of the cubic oscillator, and this angle is 
therefore chosen for all numerical calculations reported henceforth in the current investigation. The coefficients of 
the strong-coupling expansion are given in Table II, and a comparison of the values obtained by method I is made in 
Fig. 3. 

We conclude that the numerical approach (method I) can be verified to high accuracy against both weak and strong- 
coupling expansions and is found to be the most convenient method to cover all ranges of the coupling constant g. 
Furthermore, the mutual agreement of this numerical method with both the strong and the weak-coupling expansions 
confirms that an estimation of the numerical uncertainty of the results based on the apparent convergence of the 
energy levels under an increase of the size of the basis of states appears to be reliable. This observation means that 
we are, in principle, in a position to use the numerically determined resonance energies and resonance for an adiabatic 
time propagation algorithm, which will be the subject of the next section of this article. 
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III. UNIFYING STRUCTURE AND DYNAMICS 

We now attempt to reconcile the structure of the cubic potential with its dynamics in a unified framework, inspired 
by a number of investigations regarding quantum dynamics formulated in complex coordinates [32-35] . Notice that 
the construction of a time propagation algorithm from complex coordinates was explicitly mentioned as a desirable 
goal a rather long time ago, in Ref. [32]. Note also that the time propagation of wave packets in a cubic potential 
is not completely trivial: We consider a particle initially at rest and located at q(t = 0) < — (Sy/g)^ 1 in the cubic 
potential V(q) = \ q 2 + y/gq 3 - According to classical mechanics, this particle reaches q = — oo in a finite time (as is 
well known), and for quantum mechanics, this means that the component of a wave packet located, loosely speaking, 
to the left of the cubic well is accelerated toward q = — oo by the cubic term, consistent with the Ehrcnfcst theorem, 
and escapes any (necessarily finite) grid in coordinate space used for the time propagation in a finite time. This 
"escape mechanism," which leads to a loss of probability amplitude for any part of the wave packet located in a finite 
subinterval of coordinate space, affords a physically intuitive explanation for the finite decay width associated with 
the resonance eigenstates of the cubic potential. 

We proceed as follows. Using a basis spanned by the standard harmonic oscillator wavefunctions {</>j(<7)}j=0' we 
expand the cubic eigenfunctions $>n(<1) of the complex scaled Hamiltonian (13), which fulfill H c {9) $jv(<z) = En $j\r(g) 
with complex En, as follows: 

oo 

M«) = 5^cjv,j0j(g). (18) 
j=o 

Here, the complex coefficients cn.j are found by the diagonalization of the complex scaled Hamiltonian (13). This 
immediately implies that the eigenfunctions of the (dilationally-transformcd) cubic potential arc also complex. In 
Fig. 4, for example, we display the real and imaginary parts of the eigenfunctions of the ground (N = 0) and excited 
(N = 2) states of the cubic potential with ^fg = 0.1. 

After the evaluation of the eigenfunctions (18), the question may arise whether these function form a complete basis 
set. In fact, such a question is not completely trivial: it is known [35-37] that for very special values of the dilational 
parameter 9, an incomplete basis of the Hamiltonian can be obtained (sec Sec. 2.5 of Ref. [35]). In particular, it 
can be shown that for special, isolated values of 8, the number of linearly independent vectors of the complex-scaled 
potential may be smaller than the size of the Hamiltonian matrix. We note, however, that any infinitesimally small 
variation of 6 turns the spectrum into a complete one [35]. Below, therefore, we assume that the set of eigenfunctions 
(18) form a complete orthonormal basis in the complex-transformed space. That is, they fulfill the condition: 

($jv|$m) = 5 N m, (19) 

where (■]•) is the so-called c-inner product introduced by N. Moiseyev and coworkers [35, 36]: 

($jv|$ m ) = f $ N (q)$M(q)dq. (20) 

In order to illustrate the properties of the c-inner product, we represent the Schrodinger equation as a matrix eigenvalue 
problem [35] with a Hamiltonian matrix H , which wc denote in boldface notation in order to emphasize that it refers 
to the matrix obtained after the complex rotation. In such a representation, the "bra" and "ket" eigenstates of the 
Hamiltonian are given by the &f and row and column eigenvectors which satisfy the matrix equations 

H&f = Ej9ff , (21) 

and 

(Sf )* H = (if 4 #f )* = Ej (Sf )* . (22) 

We recall that the eigenvalues of a matrix and its transpose are necessarily the same (as follows from the secular 
equation) , while the corresponding left and right eigenvectors are not necessarily the same. Since the left eigenvectors 
3>j are to ipso transposed and satisfy the eigenvalue equation of the transposed Hamiltonian, it is not necessary to 
invoke complex conjugation in the definition of the inner product (20). When the Hamiltonian matrix is derived from 
an originally "purely real" Hamiltonian such as (1) by complex scaling and thus symmetric (H = H t ), then the right 
and left vectors are identical, <E> N = = *&n- This is the case for the cubic Hamiltonian. 

Wc now return to the eigenfunctions (18) of the cubic Hamiltonian. These wavefunctions, which form a complete 
basis, may be utilized for the numerical integration of the (single particle) Schrodinger equation and, hence, for 
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FIG. 4: (color online.) Real (top panels) and imaginary (bottom panels) parts of 
the wavefunction $jv of the cubic Hamiltonian as a function of coordinate q. Results 
have been obtained for the ground state N = and for the second excited state with 
N — 2. The coupling used is ^fg = 0.1 and the universal complex rotation angle is 
9 = 36°. Note that both the real and the imaginary parts of the ground-state wave 
function have nodes, but the wave function has no zero when considered as a complex 
variable. 



studying the dynamical evolution of the wave packet in the cubic potential. Note, however, that this packet is defined 
initially in the normal (i.e., not in the scaled) space. The time propagation, therefore, requires first a complex scaling 
of the initial wave packet: 

tt(g,t = 0) -> V(qc ie ,t = 0) = * c (g,t = 0) . (23) 
After the complex scaling, the initial wave packet can be expanded in the basis of the functions (18) 



# c ( ? ,f = 0)= J2 b N®N(q), 



(24) 



where = ($jv|^ c (f = 0)). The propagation of the wave packet ^f c (q,t) in the complex scaled coordinates is finally 
given by: 



(25) 



After performing the time propagation, the final function ^ c (q,t) should be transformed back into the "normal" 
space, to obtain the wave function ^(q, t). The back-transformation proceeds by a simple inversion of the replacement 
operation (23), 



y c (q,t)^V(q,t) = V c (qe- i6 ,t) 



(26) 



Using Eq. (25), we perform calculations for the time propagation of an (initially) Gaussian wave packet within the 
cubic potential with the coupling parameter yfg = 0.04. In order to visualize the results of such a time propagation, 
we display in Fig. 5 the so-called autocorrelation function 



(27) 
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Time 



FIG. 5: (color online.) Time evolution of the autocorrelation function (27). Calcula- 
tions are performed for an (initially) Gaussian wave packet of the form *&(£/, f = 0) = 
e~ q /(71-/2) 1 / 4 in a cubic potential with the coupling constant ^fg = 0.04. Results 
obtained from Eq. (25) are compared with those of the well-known Crank-Nicolson 
approach (solid line). The lower panel is a close-up of the upper one, where the range 
of the ordinate axis in restricted to the interval P(t) £ (0.98, 1.0). The maximum of 
P(t) — 0.9939 near the sixth oscillation, at time t = 19.29, is faithfully reproduced 
by the complex scaled propagation. 



which is commonplace for studying tunneling phenomena in multi-well quantum dot potentials [38]. As seen from 
Fig. 5, the wave packet, initially located inside the well (at q = 0), performs oscillations confined within the well. 
With every cycle of the oscillation, however, the autocorrelation function slightly decreases due to tunneling. We have 
also verified numerically that the decay of a wave packet corresponding to a single resonance eigenstate is described 
correctly by the back-transformation (26), via a comparison to a Crank-Nicolson method, where the spatial box had to 
be chosen sufficiently large in the latter case to accommodate the spreading of the wave packet for larger propagation 
times. 

We note that suitable generalizations of the method (25) can also be applied to time-dependent (driven) potentials 
with resonances, e.g., within the adiabatic approximation, which is valid provided one uses sufficiently small time steps 
so that the potential can be regarded as constant in time within each time step. Recently, the adiabatic approach has 
been successfully applied to study the propagation of the wave packet in the driven double-well-like potential [20] . 

Up to now, we have considered the idealized case of a decoherence-free, dissipation-free time evolution of the 
metastable states in potentials which can be approximated by a cubic Hamiltonian. The coupling of the system to 
an ensemble of harmonic oscillator modes, whose levels are occupied according to thermal distributions, changes the 
dynamics. The latter model has been studied extensively in the literature where it is commonly referred to as the 
Caldcira-Lcggett model [39] , and it has been discussed at length in a number of research articles and books [40-48] ; 
therefore its details are relegated to the Appendix A. The model describes dissipation and the corresponding quantum 
decoherence, and, in general, leads to multiplicative corrections to the decay widths of the resonance eigenstates, of 
the form 

F N -> T N (1 + S N ) (28) 

for the decay widths. The corrections Sjy, initially, could be assumed as affecting only the structure of the resonance 
energies, but using our time propagation, it is easy to include them into the dynamics as well, by simply applying the 
replacement (28) to the Tn which are present in Eq. (25). 



IV. SUMMARY AND OUTLOOK 



In this paper, we have illustrated a certain duality of the cubic anharmonic oscillator; for purely imaginary coupling, 
the eigenenergies are real, whereas for real coupling, the resonance energies are complex (see Sees. II A and II B). 
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Moreover, for imaginary coupling, there are no instanton configurations to consider, and the quantization condition 
involves only the "perturbative B function," suitably resummed, and is of the plain Bohr-Sommerfeld type [see 
Eq. (7)]. By contrast, for real coupling, there are instanton configurations present, and these manifest themselves 
in nonvanishing decay widths of the state and in a modified quantization condition with allowance for the instanton 
configurations. The modified quantization is not obtained here, but we lay the groundwork for its construction. 
[Specifically, the leading term of the "instanton A function" is given in Eq. (12).] 

In other words, even anharmonic oscillators (quartic, sextic, etc.), for real coupling, arc described by a plain Bohr- 
Sommcrfeld quantization condition, whereas for negative coupling, they give rise to instantons. The structure of an 
even oscillator for negative coupling is thus in a certain sense analogous to an odd oscillator for real coupling ^fg £ R, 
and the presence of instantons and the corresponding nontrivial saddle points of the Euclidean action demand a 
modification of the Bohr-Sommerfeld quantization condition. By contrast, an even oscillator for positive coupling 
parameter is analogous to an odd oscillator for purely imaginary coupling ^fg = i(3 and is described by a plain 
Bohr-Sommerfeld quantization condition and completely characterized by Rayleigh-Schrodinger perturbation theory. 

For the determination of the resonance energies of the cubic potential (positive coupling, see Sec. II B), three different 
methods have been compared. Method I (numerical approach in complex coordinates) is found to be more universal 
in applicability than method II (resummed weak-coupling expansion) and method III (strong-coupling expansion), 
although good agreement is observed in the specified domains of applicability for each of the latter two methods. In 
any case, the numerical approach provides us with resonance energies and corresponding wave functions in complex 
space which can be used in order to unify the analysis of the structure and of the dynamics of quantal particles. 
Hereby, the immediate physical interpretation of the decay widths can be illustrated in a particularly clear way by 
investigating the time propagation of wave packets. 

As shown in Sec. Ill, the complex rotated eigenfunctions evolve in time according to complex resonance energies 
which describe the quantum tunneling and the decay of the wave packets. When transformed back into real space, 
the algorithm leads to results which are in agreement with traditional quantum dynamics solvers like the Crank- 
Nicolson method. The structure and the dynamics of the potentials are thus obtained in a single, unified method 
which allows for an inclusion of modifications to the widths of the resonances due to dissipation (coupling to the 
environment). For large propagation times, our algorithm is numerically stable, because the time evolution of the 
resonance eigenstates given by Eq. (25) remains valid for the entire coordinate space. The "escape" of the wave packet 
toward q — > — oo, which is more pronounced for the cubic potential than, e.g., for the linear ponderomotive potential 
known from laser physics, is thus automatically incorporated into our algorithm. Moreover, slight modifications of the 
structural properties of the resonances due to additional perturbations (e.g., couplings to an environment) can easily 
be incorporated into the dynamics, as exemplified by Eq. (28) and discussed in more detail in Appendix A below. 
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APPENDIX A: DISSIPATION IN QUANTUM TUNNELING 

Recently, a rather large number of theoretical [49-52] as well as experimental [53-57] studies have been devoted 
to the current-biased Josephson junction since its low-lying energy levels can be used to implement a solid-state 
quantum bit (qubit). This qubit circuit might be considered as a conceivable basis of a scalable quantum computer. 
The theoretical analysis of the structure and dynamics of Josephson junction qubits can be traced back to the cubic 
potential, at least in the limit of a dissipation-free, or decoherence-free environment, which is an idealized scenario 
that is nevertheless pursued in the construction of quantum computing devices. The measurement of the qubit state 
utilizes the escape from the cubic potential via tunneling and thus it is important to have the dynamical aspects of 
related potential under control. 

When we work with a plain cubic Hamiltonian, it is clear that the analysis applies only to the idealized case of a 
decoherence-free, dissipation-free time evolution of metastable states in potentials which can be approximated by a 
cubic Hamiltonian. The coupling of the system to an ensemble of harmonic oscillator modes, whose levels arc occupied 
according to thermal distributions, changes the dynamics and can be incorporated approximately into our algorithm 
by the replacement (28). The so-called Caldeira-Leggett model has been studied extensively in the literature [39-48] 
and describes dissipation and the corresponding quantum decoherence. Ideally, of course, one has no dissipation in 
a device used for quantum computing, and this scenario has been the focus of the current investigation. It might 
be worth emphasizing that our approach allows for a full access to the quantum dynamics in the limit of vanishing 
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dissipative terms, whereas the approach originally outlined in Rcfs. [39-43] can be used to investigate the rate at 
which the system loses information about coherent quantum states in the presence of dissipation, but does not lead 
to a general description of the full quantum dynamics in timc-dcpcndcnt potentials. 

In general, the coupling to other degrees of freedom leads to the phenomenon of dissipation where the energy is 
transferred irreversibly from the system under consideration to the environment. In the work by A. O. Caldeira and 
A. J. Leggett [39], the so-called "system-plus-bath" model has been introduced to describe dissipation. Within this 
model, the bath (environment) is considered to be representable as a set of harmonic oscillators interacting linearly 
with the system under consideration (i.e., with our cubic potential). The total Hamiltonian for such a system can be 



written as 



^ / 2 2 2 \ 

+ W^ + *^ + ^ * \ (A1) 

where the first line displays "environmentally decoupled" system, i.e. the cubic anharmonic oscillator in our case. The 
sum over i in Eq. (Al) contains the Hamiltonians for a set of N harmonic oscillators which arc bilincarly coupled 
with strength Cj to the system. Finally, the last term represents a potential renormalization term [47, 48] . 

For the practical application of Eq. (Al), it is necessary to eliminate the external degrees of freedom, i.e. the bath 
modes labeled 1, . . . , N. Since these are just harmonic oscillator modes, we can easily solve the equations of motion 
for the external degrees of freedom. In a classical framework, we can derive the equation of motion for the cubic 
anharmonic oscillator in the bath as [47, 48] 

/ 

2 



q + I dt'i(t - t')q(t') + q + 3^gq 2 = £(t) . (A2) 
Here, the damping kernel 7 is given by 



N 



■2 



= Y, cosM) , (A3) 



i=i 



and the fluctuating acceleration 

m = -7(*)«(o) 



N 

E 



d I Xi(0) cos(wit) + sin(cJii) I (A4) 



m l uj l 



i=i 

contains the initial conditions for the bath modes. 

On a fully quantum level, it is convenient to describe the decay of the metastable state coupled to environment 
within the path integral formalism, because harmonic oscillator modes can easily be integrated out. The analysis 
of dissipative quantum tunneling is naturally performed in the Euclidean space (see Refs. [47, 48] for more details). 
That is, after a suitable Wick rotation, the imaginary-time path integral provides a representation of the equilibrium 
density matrix which can be represented as 

Pf>{Qi,Qf) = Y- J V q c^(-S E [q]) , (A5) 

g(o)=«i 

where Zp is the partition function which here acts as a normalization prefactor, and (3 = l/iksT) is the imaginary or 
thermal time (we set % = 1). The paths in the integral (A5) are weighted with a phase factor that contains the effective 
Euclidean action for the anharmonic oscillator obtained after integrating out the modes of thermal bath [39, 40, 47, 48] 
and reads 



P 



S E [q] = JdT^ + ^ + ^q 3 


P 

+ I I dr [ dr'k(T - t') q(r) q(T') . (A6) 



2 
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The first term of this action is just a standard Euclidean action of the cubic Hamiltonian (1), while the second term 
describes the frictional influence of the environment. The temperature-dependent, nonlocal (in time) damping kernel 
k(r — t') may be expressed as a Fourier series (see, e.g., Refs. [44, 48]) 

+oo 

k{ f) = pYl &»exp(ii/„t), (A7) 

n=— co 

where the v n = 27m//3 are bosonic Matsubara frequencies. Finally, we have £ n = \v n \ 7(|^„|), where 7 denotes the 
Laplace transform of the damping kernel (A3). 

Based on the path integral approach (A5)-(A7), the quantum-mechanical decay of the cubic oscillator under the 
influence of a heat bath environment has been studied. In particular, explicit results have been obtained which helped 
to understand how the tunneling rate depends on the temperature of the environment, in the presence of medium 
and strong damping terms [41, 44, 48]. 
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